rm(list=ls())
#setwd("/Main Text-Figures/Figure1/")

#install.packages("pscl")
#install.packages("emIRT")
library(pscl)
library(emIRT)

########################################

#simulation to create:
  
  ##"idealpoints-fig1a-data.RData" k <- .3
  ##"idealpoints-fig1b-data.RData" k <- .2
  ##"idealpoints-fig1c-data.RData" k <- .1

alpha <- .5
k <- 0.1 #0.2, 0.3
N <- 100
v <- 1000
nsim <- 25

  phi <- runif(1,.2,.8)
  
  b1 <- runif(v,0,1)
  b2 <- b1
  sq1 <- runif(v,0,1)
  sq2 <- sq1
  
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  
  d1 <- runif(1,.5,2)
  d2 <- d1
  
########################################
  
  hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
  colnames(hypo.pro.1d) <- c("dim","d","b","sq")
  hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
  colnames(hypo.pro.2d) <- c("dim","d","b","sq")
  hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
  hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
  
########################################
  
  xL1 <- c(0,runif(npL-1,min=0-d1,max=d1))
  xL2 <- xL1
  xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) 
  xR2 <- xR1
  
  left.party <- as.data.frame(cbind(round(xL1,3),round(xL2,3),rep(0,npL)))
  colnames(left.party) <- c("x1","x2","partyid")
  right.party <- as.data.frame(cbind(round(xR1,3),round(xR2,3),rep(1,npR)))
  colnames(right.party) <- c("x1","x2","partyid")
  all.leg <- as.data.frame(rbind(left.party,right.party))
  
########################################
#votes for it?
  
  vote.func1 <- function(x1,b1,sq1){
    as.numeric(abs(x1-b1)  < abs(x1-sq1) )
  }
  vote.func2 <- function(x2,b2,sq2){
    as.numeric(abs(x2-b2)  < abs(x2-sq2) )
  }
  out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
  for(i in 1:N){
    for(j in 1:v){
      out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
    }
  }
  for(i in 1:N){
    for(j in (v+1):(2*v)){
      out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
    }
  }
  
########################################
#pass?
  
  passes <- as.numeric(colSums(out) > .5*N)
  hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
  for (i in 1:nrow(hypo.pro.all)){
    hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
  }
  
########################################
#legislator unity scores

  left.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
  }
  right.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
  }
  out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
  
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] - right.aligned[j]
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]-left.aligned[j]
      } else {
        0
      }
    }
  }
  
  out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] 
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]
      } else {
        0
      }
    }
  }
  
  ########################################
  #roll call votes requested?
  
  out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
  for(i in 1:N){
    for(j in 1:(v*2)){
      out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
    }
  }
  hypo.pro.all$left.cscore <- out.cohesion[1,] 
  hypo.pro.all$right.cscore <- out.cohesion[N,]
  
  hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
  hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
  
  hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
  hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
  
  hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)

########################################
#proposed?
  
  for(i in 1:nrow(hypo.pro.all)){
    hypo.pro.all$proposing.party[i] <- if(hypo.pro.all$b[i] < hypo.pro.all$sq[i]) 0 else 1 
  }
  
  hypo.pro.all$proposed.left <- rep(0,length(nrow(hypo.pro.all)))
  for(i in 1:nrow(hypo.pro.all)){
    if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else {
      0
    }
  }
  
  hypo.pro.all$proposed.right <- rep(0,length(nrow(hypo.pro.all)))
  for(i in 1:nrow(hypo.pro.all)){
    if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else {
      0
    }
  }
  
  hypo.pro.all$proposed <- hypo.pro.all$proposed.left+hypo.pro.all$proposed.right
  all.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1),]
  rcv.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1 & hypo.pro.all$rcv.request==1),]

########################################
  #vote matrixes of proposed votes
  vote.matrix.all <- out[,all.proposals$id]
  vote.matrix.rcv <- out[,rcv.proposals$id]

########################################
  
  y.all <- as.data.frame(vote.matrix.all)
  rc.all <- rollcall(y.all)
  rc.all$votes[rc.all$votes==0]<- -1
  rc.all$codes$nay <- -1
  rcfast.all <- convertRC(rc.all)
  p.all <- makePriors(rcfast.all$n, rcfast.all$m, 1)
  s.all <- getStarts(rcfast.all$n, rcfast.all$m, 1)
  
  fast.out.all <- binIRT(.rc = rcfast.all,
                         .starts = s.all,
                         .priors = p.all,
                         .control = {
                           list(threads = 6,
                                verbose = FALSE,
                                thresh = 1e-6
                           )
                         },
                         .anchor_subject=100
  )
  full.est.all <- round(fast.out.all$means$x,3)
  full.est.all.party <- cbind(seq(1:length(fast.out.all$means$x)),full.est.all,all.leg)
  colnames(full.est.all.party) <- c("id","full.est.irt","idp1.true","idp2.true","party.id")
  head(full.est.all.party)
  
  y.rcv <- as.data.frame(vote.matrix.rcv)
  rc.rcv <- rollcall(y.rcv)
  rc.rcv$votes[rc.rcv$votes==0]<- -1
  rc.rcv$codes$nay <- -1
  rcfast.rcv <- convertRC(rc.rcv)
  p.rcv <- makePriors(rcfast.rcv$n, rcfast.rcv$m, 1)
  s.rcv <- getStarts(rcfast.rcv$n, rcfast.rcv$m, 1)
  
  fast.out.rcv <- binIRT(.rc = rcfast.rcv,
                         .starts = s.rcv,
                         .priors = p.rcv,
                         .control = {
                           list(threads = 6,
                                verbose = FALSE,
                                thresh = 1e-6
                           )
                         },
                         .anchor_subject=100
  )
  full.est.rcv <- round(fast.out.rcv$means$x,3)
  full.est.rcv.party <- cbind(seq(1:length(fast.out.rcv$means$x)),full.est.rcv,all.leg)
  colnames(full.est.rcv.party) <- c("id","rcv.est.irt","idp1.true","idp2.true","party.id")
  head(full.est.rcv.party)
  
  dat <- full.est.rcv.party
  dat$full.est.irt <- full.est.all.party$full.est.irt
  colnames(dat) <- c("id","partial","true","true2","party","full")
  
########################################
  
  dat.order <- dat[order(dat$true),]
  dat.order <- as.data.frame(cbind(seq(1:N),dat.order))
    colnames(dat.order) <- c("true.rank","id","partial","true","true2","party","full")
  
  dat.left <- dat.order[which(dat.order$party==0),]
  dat.right <- dat.order[which(dat.order$party==1),]
  
########################################

#figures 1a, 1b, 1c
load("idealpoints-fig1a-data.RData")
load("idealpoints-fig1b-data.RData")
load("idealpoints-fig1c-data.RData")
  
#pdf("Fig1c_.pdf",width=6,height=5)
  par(mar=c(4,4,4,1))
  plot(dat.left$true,dat.left$full,pch=15,col="grey80",ylim=c(-4,3),xlab="",ylab="",xaxt='n',yaxt='n',cex=1.75,bg="white",xlim=c(-2,3)) 
  points(dat.right$true,dat.right$full,pch=16,col="black",cex=1.75,bg="White")
  points(dat.left$true,dat.left$partial,col="grey60",pch=22,bg="white",cex=1.25) 
  points(dat.right$true,dat.right$partial,col="black",pch=21,bg="white",cex=1.25) 
  mtext(side=1,line=2.5,"True Ideal Points",cex=1.5)
  mtext(side=2,line=2.5,"Estimated Ideal Points",cex=1.5)
  axis(side=1,at=seq(-2,3,.5),labels=c("-2","","-1","","0","","1","","2","","3"))
  axis(side=2,at=seq(-4,3,1),labels=c("-4","","-2","","0","","2",""),las=1)
  mtext(side=3,line=2,print(paste("RCV on",ncol(y.rcv),"of",ncol(y.all),"Votes")),cex=1.5)
  mtext(side=3,line=1,print(paste("d=",round(d1,2),"; phi=",round(phi,2),"; k=",k)))
#dev.off()
  
#pdf("legend.pdf",width=14,height=.5)
#  par(mar=c(0,0,0,0))
#  plot(0,0,xlim=c(-1,1),ylim=c(.5,1),col="white",ylab="",xlab="",xaxt='n',yaxt='n',bty='n')
#  legend(-1,1,legend=c("Left Party, Full Sample","Left Party, RCV Sample","Right Party, Full Sample",
#                                            "Right Party, RCV Sample"),pch=c(15,22,16,21),
#                               col=c("grey60","grey60","black","black"),cex=1.5,bty='n',bg='n',ncol=4)     
#dev.off()